function [rotation, Fh_torque,Fhs] = adaptiveLoadPath(data,opts,plotThings)




%%% ------- set default values ---------------------------
dFh0 = 2; %heel pull force increment
exitTol = 1e-1; %exit value for current slope of (dF/da) over initial dF/da
maxStepIts = 5000; %max number of steps to follow load path
maxTryIts = 100; %max number of iterations to decrease load increment to get 
% around non-linearities and horriblness
maxToeAngle = 45*pi/180;

a0 = [0.;0]; % initial toe and bellows angle
options.a = 0; %set dummy options set
options.maxIts = 1e2; %max number of newton iterations
relax0 = 1.0; %default relaxation
relax2 = 0.9; %relaxation to use when bad point found and step size decreasess
options.relax =relax0; %relaxation for newton solver, should be between

%%% set values if not input, check if options were input
if nargin < 3
    plotThings = 0;
end

if nargin > 1
    if isfield(opts,'dfh')
        dFh0 = opts.dfh;
    end
end



%%% actually do things
if plotThings
    hh = figure;
    axis equal
    grid on
    disp('adjust boot visualization plot for viewing during sim, then hit enter')
    pause
end

tic

nVals0 = 1e3; %number of values to intialize storage arrays with

Fh_torque = zeros(1,nVals0);
rotation = zeros(1,nVals0);
aNs = zeros(2,nVals0);
Fhs = zeros(1,nVals0);


%figure out initial heel lift force 
[m0] = basicBoot(a0,0,data);
Fh0 = m0(1)/data.boot.soleLength;
[alphaN,sigM] = dampedNewton(@basicBoot,a0,options,abs(Fh0),data);
[~,~,~,Fh_torque(1),rotation(1)] = basicBoot(alphaN,abs(Fh0),data);
aNs(:,1) = alphaN;
Fhs(1) = abs(Fh0);

dFh = dFh0;
it = 1;
done = 0;
while ~done
    it = it+1;
    
    succeeded = 0;
    tryIt = 0;
    while ~succeeded
        tryIt = tryIt + 1;
        
        Fhs(it) = Fhs(it-1) + dFh;
        
        [alphaN,sigM,~,dumped] = dampedNewton(@basicBoot,alphaN,options,Fhs(it),data);
        
        if dumped
            dFh = dFh/2;
            options.relax = relax2;
            disp('Dropped dFh')
        else
            options.relax =relax0;
            dFh =dFh0;
            succeeded =1;
        end
        
        if tryIt > maxTryIts
            error('Exceeded adaption iterations')
        end
        
    end
            

    
     [~,~,~,Fh_torque(it),rotation(it)] = basicBoot(alphaN,Fhs(it),data);
     
     aNs(:,it) = alphaN;
     
     if any(isnan(alphaN))
         Fh_torque = Fh_torque(1:end-1);
         rotation =rotation(1:end-1);
         aNs = aNs(:,1:end-1);
         break;
     end
     
     if it ==2
         initialSlope = max(abs((aNs(:,2) - aNs(:,1))./(Fhs(2) - Fhs(1))));
     elseif it > 3
         cSlope = max(abs((aNs(:,it) - aNs(:,it-1))./(Fhs(it) - Fhs(it-1))));
         
%          abs(cSlope/initialSlope)
         
         if abs(cSlope/initialSlope) < exitTol
             done = 1;
         elseif maxStepIts < it
             error('Exceeded number of force incremenet iterations')
         elseif alphaN(1) > maxToeAngle
             done = 1;
         end
     end
         
     
     
     if plotThings
         clf
        axis equal
        grid on
         plotBindingPosition(data,alphaN,hh,'b')
         tstring = ['[FhL,alpha-toe (deg),alpha-bellows (deg)]=[',num2str(Fhs(it),'%2.3g'),' , ',...
             num2str(alphaN(1)*180/pi,'%2.3g'),' , ',num2str(alphaN(2)*180/pi,'%2.3g'),']'];
         title(tstring)
     
     
        pause(.01)
        
     end
end

rotation = rotation(1:it);
Fh_torque = Fh_torque(1:it);
aNs =aNs(:,1:it);
Fhs = Fhs(1:it);

toc

rotation = rotation*180/pi;

if plotThings
    figure
    plot(rotation,Fh_torque,'.-')
    xlabel('Heel rotation angle')
    ylabel('Heel lift torque')
    grid on
end